program afgeleide
	use ieee_exceptions    
	implicit none
    real, parameter :: ALPHA = 1
    integer, parameter :: wp = kind(0d0)
    real(wp) :: startwaarde = 1
    real(wp) :: approx, exact, relerr
    integer :: i
    logical :: flag_value
    call ieee_set_halting_mode(ieee_divide_by_zero, .false.)
    do i=1,16
        call ieee_set_flag(ieee_divide_by_zero, .false.)
        approx = approx_df(startwaarde,10.0_wp**(-i))
        exact = df(startwaarde)
	relerr = (approx-exact)/exact
    	call ieee_get_flag(ieee_divide_by_zero, flag_value)
        if(flag_value) then
        	write(*, fmt='("h=10^",i3,", approx=",es9.2,", exact=",es9.2,"-> Divide by zero")'),-i,approx,exact
		relerr = 1.0
	else
        	write(*, fmt='("h=10^",i3,", approx=",es9.2,", exact=",es9.2,"-> 			relerr=",es9.2)'),-i,approx,exact,relerr
	endif
	call ieee_set_flag(ieee_divide_by_zero, .false.)
    enddo

contains

function f(x) result(value)
    real(wp), intent(in) :: x
    real(wp) :: value
    value = (x-ALPHA)**(-7)
end function

function approx_df(x,h) result(value)
    real(wp), intent(in) :: x, h
    real(wp) :: value
    logical :: overflow, underflow, inexact
    call ieee_set_flag(ieee_overflow, .false.)
    call ieee_set_flag(ieee_underflow, .false.)
    call ieee_set_flag(ieee_inexact, .false.)
    value = (f(x+h)-f(x-h))/(2*h)
    call ieee_get_flag(ieee_inexact, inexact)
    if(inexact) print *, 'inexact'
    call ieee_get_flag(ieee_overflow, overflow)
    if(overflow) print *, 'overflow'
    call ieee_get_flag(ieee_underflow, underflow)
    if(underflow) print *, 'underflow'
    call ieee_set_flag(ieee_overflow, .false.)
    call ieee_set_flag(ieee_underflow, .false.)
    call ieee_set_flag(ieee_inexact, .false.)
end function

function df(x) result(value)
    real(wp),intent(in) :: x
    real(wp) :: value
    value = -7*(x-ALPHA)**(-8)
end function

end program
    
